%平台集成
function [location_detect,location_real,LocationError_avg,T,point1,point2]=fun_2(filepath2_DATA,filepath2_HEX,filepath_locationReal,Area)
T=0;
%% 读取仿真后数据
data1=load(filepath2_DATA);
data2=load(filepath2_HEX);
data3=load(filepath_locationReal);
DATA = data1.DATA;
HEX = data2.HEX;
location_real=data3.location_real;
%% 数据区域选择
if Area==1
    lat1=30;lat2=34;lon1=-99.5;lon2=-95.5;
end
if Area==2
    lat1=42;lat2=46;lon1=-116;lon2=-112;
end
if Area==3
    lat1=28;lat2=32;lon1=-101;lon2=-97;
end 
%% 存储每条轨迹的对应数据，筛选出数据点超过100的数据
H=cell(500,4); %元胞H存储每条轨迹的对应信息
uni = unique(HEX); %筛选不同的HEX
k=1;
for i = 1:1:length(uni)
    u = uni(i); %u是每一个HEX号
    id = find(strcmp(HEX,u)==1);    %找到数组HEX[]中与u相同的元素的位置,存入数组id[]中
    if length(id)>=100 
        H{k,1}=u; %储存航班号
        H{k,2}=id; %储存在原数据中的位置
        H{k,3}=DATA(id,6); %储存时间信息
        H{k,4}=[DATA(id,1),DATA(id,2),DATA(id,7),DATA(id,8),DATA(id,3)]; %储存位置信息:经度、纬度、高度
        H{k,5}=DATA(id,4);%储存速度信息
        k=k+1;
    end
end
H=H(1:k-1,:); % 轨迹数量k-1

%% 找到发生断裂的轨迹track_hex,找到消失和出现的位置
time=unique(DATA(:,6)); %time为不重复的时间序列
track_hex=strings; %track_hex是一个字符串数组，储存异常航班号
track_locat=zeros(1,13);
k=1;
 for i = 1:1:length(H) %遍历每一条轨迹
    for j=1:1:length(H{i,2})-1 %轨迹上的一个点j
        for a=1:1:length(time)-1 %找到当前时间节点a
            if H{i,3}(j,1)==time(a)
               break;
            end
        end
        for a_1=a+1:1:length(time) %找到下一个时间节点a_1
            if H{i,3}(j+1,1)==time(a_1)
                break;
            end
        end
        if a_1-a>30 && a_1-a<360%即断裂时间>5分钟 && <1小时
            if sqrt((H{i,4}(j,3)-H{i,4}(j+1,3))^2+(H{i,4}(j,4)-H{i,4}(j+1,4))^2)>20 %两点之间距离>0.6
            track_hex(k,1)=H{i,1}; %track_hex存储断裂hex:0d0efa
            track_locat(k,1)=i; %储存hex的编号:2
            track_locat(k,2)=H{i,4}(j,1);  %储存轨迹断裂点的经度
            track_locat(k,3)=H{i,4}(j,2);  %储存异常断裂点的纬度
            track_locat(k,4)=H{i,4}(j,3);  %储存异常断裂点的横坐标
            track_locat(k,5)=H{i,4}(j,4);  %储存异常断裂点的纵坐标
            track_locat(k,6)=H{i,4}(j,5);  %储存异常断裂点的高度
            track_locat(k,7)=H{i,4}(j+1,1);  %储存轨迹恢复点的经度
            track_locat(k,8)=H{i,4}(j+1,2);  %储存轨迹恢复点的纬度
            track_locat(k,9)=H{i,4}(j+1,3);  %储存轨迹恢复点的横坐标
            track_locat(k,10)=H{i,4}(j+1,4);  %储存轨迹恢复点的纵坐标
            track_locat(k,11)=H{i,4}(j+1,5);  %储存异常断裂点的纵坐标
            k=k+1;
            end
        end
    end
 end

if isempty(track_locat) %如果为空
    T=1;
    location_detect=[];LocationError_avg=[];point1=[];point2=[];
    return 
end

%% 高度<200
i=1;
while i<=size(track_locat,1)
    if track_locat(i,6)<200
        track_locat(i,:)=[];
        track_hex(i,:)=[];
        i=i-1;
    end
    i=i+1;
end
if isempty(track_locat) %如果为空
    T=2;
    location_detect=[];LocationError_avg=[];point1=[];point2=[];
    return 
end

 %% DBSCAN聚类
%断裂点聚类
eps=30; %%规定聚类半径
MinPts=8; %%规定半径内最小数目点
point1=[];point1(:,1)=track_locat(:,4);point1(:,2)=track_locat(:,5);point1(:,3)=track_locat(:,6);
IDX_1=DBSCAN(point1(:,1:2),eps,MinPts);%调用DBSCAN函数，输入参数有三个：原始数据，聚类半径，最小数目点
track_locat(:,12)=IDX_1;

%% 恢复点聚类
eps=30;%%规定聚类半径
MinPts=8;%%规定半径内最小数目点
point2=[];point2(:,1)=track_locat(:,9);point2(:,2)=track_locat(:,10);point2(:,3)=track_locat(:,11);
IDX_2=DBSCAN(point2(:,1:2),eps,MinPts);%调用DBSCAN函数，输入参数有三个：原始数据，聚类半径，最小数目点
track_locat(:,13)=IDX_2;

%% 排除噪声点，两种聚类结果中的噪声点
i=1;
while i<=size(track_locat,1)
    if track_locat(i,12)==0 || track_locat(i,13)==0
        track_locat(i,:)=[];
        track_hex(i,:)=[];
        IDX_1(i)=[];
        IDX_2(i)=[];
        i=i-1;
    end
    i=i+1;
end
if isempty(track_locat) %如果为空
    T=3;
    location_detect=[];LocationError_avg=[];point1=[];point2=[];
    return
end
%
track_locat(:,2)=track_locat(:,4);
track_locat(:,3)=track_locat(:,5);
track_locat(:,4)=track_locat(:,6);
track_locat(:,5)=track_locat(:,9);
track_locat(:,6)=track_locat(:,10);
track_locat(:,7)=track_locat(:,11);

track_locat(:,8)=[];
track_locat(:,8)=[];
track_locat(:,8)=[];
track_locat(:,8)=[];

%
point1=[];point2=[];
point1(:,1)=track_locat(:,2); %距离坐标下的断裂点和恢复点
point1(:,2)=track_locat(:,3);
point1(:,3)=track_locat(:,4);
point2(:,1)=track_locat(:,5);
point2(:,2)=track_locat(:,6);
point2(:,3)=track_locat(:,7);
%% 分出断裂点的各类点集
%分point1
%分出断裂点的各类点集
clusterData1={};
for i = 1:max(IDX_1)
    clusterData1{i,1}=track_locat(IDX_1 == i,:); %find()函数查找与条件相匹配的元素的索引
end
%% 依次验证，找到最佳分割法
clusterData1=sequence(clusterData1,1); %按断裂点排序

circle={}; circle_1={};%存储断裂点拟合圆数据
for j=1:1:size(clusterData1,1)
    circle_2={};
    data=clusterData1{j,1};
    [circle_2]=cirlce(data,1); %circle存储拟合圆数据
    t=aver_m(circle_2);  %计算误差参数m:t1
    for i=1:size(data,1)
        data=circshift(data,1); %数组后移一位
        [circle_1]=cirlce(data,1); %circle_1存储拟合圆数据
        t1=aver_m(circle_1); %计算误差参数m:t1
        if size(circle_1,1)<size(circle_2,1) %优化
            circle_2=circle_1;
            t=t1;
        end
        if size(circle_1,1)==size(circle_2,1) & t1<t
            circle_2=circle_1;
            t=t1;
        end
    end
    circle=[circle;circle_2];
end
%判断是否能合并,并储存合并后数据至circle_1
circle_1={};
[circle_1]=cb(circle,1);
% %% 
% [p,m]=Circle_Fitting(circle_1{1,1}(1:23,2:3))
% %% 
% 
% plot(circle_1{1,1}(4:12,2),circle_1{1,1}(4:12,3),'*b')
%% 依次删除个别异常点,合并计算拟合参数
%删除点数少于5的点簇
i=1;
while i<=size(circle_1,1)  
    if size(circle_1{i,1},1)<5
        circle_1(i,:)=[];
        i=i-1;
    end
    i=i+1;
end
%去头尾,计算距离参数
for t=1:1:size(circle_1,1)
    len=size(circle_1{t,1},1);
    [p,~]=Circle_Fitting(circle_1{t,1}(2:len-1,2:3));
    data=circle_1{t,1};
    r=p(3);
    a=p(1);
    b=p(2);
    for j=1:size(data,1)
        s=abs(pdist([[a,b];data(j,2:3)])-r);
        circle_1{t,1}(j,8)=s;
    end  
end
%删除阈值s>2的点
for t=1:size(circle_1,1)
    data=circle_1{t,1};
    i=1;
    while i<=size(data,1)
        if data(i,8)>2 %删除阈值s>5的点
            data(i,:)=[];
            i=i-1;
        end
        i=i+1;
    end
    circle_1{t,1}=data;
end
% 重新计算拟合参数
for t=1:1:size(circle_1,1)
    [p,m]=Circle_Fitting(circle_1{t,1}(:,2:3));
    circle_1{t,2}=p;
    circle_1{t,3}=m;
end
%判断是否能合并,并储存合并后数据至circle_1
circle_0={};
[circle_0]=cb(circle_1,1);
circle_1=circle_0;
%% 删除半径小于30的点簇,个数小于10的点簇
i=1;
while i<=size(circle_1,1)  
    if circle_1{i,2}(3)<30
        circle_1(i,:)=[];
        i=i-1;
    end
    i=i+1;
end

i=1;
while i<=size(circle_1,1)  
    if size(circle_1{i,1},1)<10
        circle_1(i,:)=[];
        i=i-1;
    end
    i=i+1;
end

%% 分出恢复点的各类点集
%分point2
%分出断裂点的各类点集
clusterData2={};
for i = 1:max(IDX_2)
    clusterData2{i,1}=track_locat(IDX_2 == i,:);% find()函数查找与条件相匹配的元素的索引
end
%% 依次验证，找到最佳分割法
clusterData2=sequence(clusterData2,2); %按恢复点排序
%
circle={};circle_3={}; %存储断裂点拟合圆数据
for j=1:1:size(clusterData2,1)
    circle_2={};
    data=clusterData2{j,1};
    [circle_2]=cirlce(data,2); %circle存储拟合圆数据
    t=aver_m(circle_2);  %计算误差参数m:t1
    for i=1:size(data,1)
        data=circshift(data,1); %数组后移一位
        [circle_3]=cirlce(data,2); %circle_3存储拟合圆数据
        t1=aver_m(circle_3); %计算误差参数m:t1
        if size(circle_3,1)<size(circle_2,1) %优化
            circle_2=circle_3;
            t=t1;
        end
        if size(circle_3,1)==size(circle_2,1) & t1<t
            circle_2=circle_3;
            t=t1;
        end
    end
    circle=[circle;circle_2];
end
%判断是否能合并,并储存合并后数据至circle_2
circle_0={};circle_2={};
[circle_2]=cb(circle,2);
%% 回代，依次删除个别异常点,合并计算拟合参数
%删除点数少于5的点簇
i=1;
while i<=size(circle_2,1)  
    if size(circle_2{i,1},1)<5
        circle_2(i,:)=[];
        i=i-1;
    end
    i=i+1;
end
%去头尾,计算距离参数
for t=1:1:size(circle_2,1)
    len=size(circle_2{t,1},1);
    [p,~]=Circle_Fitting(circle_2{t,1}(2:len-1,5:6));
    data=circle_2{t,1};
    r=p(3);
    a=p(1);
    b=p(2);
    for j=1:size(data,1)
        s=abs(pdist([[a,b];data(j,5:6)])-r);
        circle_2{t,1}(j,8)=s;
    end  
end
%删除阈值s>2的点
for t=1:size(circle_2,1)
    data=circle_2{t,1};
    i=1;
    while i<=size(data,1)
        if data(i,8)>2 %删除阈值s>5的点
            data(i,:)=[];
            i=i-1;
        end
        i=i+1;
    end
    circle_2{t,1}=data;
end
% 重新计算拟合参数
for t=1:1:size(circle_2,1)
    [p,m]=Circle_Fitting(circle_2{t,1}(:,5:6));
    circle_2{t,2}=p;
    circle_2{t,3}=m;
end
%判断是否能合并,并储存合并后数据至circle_
circle_0={};
[circle_0]=cb(circle_2,2);
circle_2=circle_0;

%% 删除半径小于30的点簇
i=1;
while i<=size(circle_2,1)  
    if circle_2{i,2}(3)<30
        circle_2(i,:)=[];
        i=i-1;
    end
    i=i+1;
end
%删除点个数小于10的点簇
i=1;
while i<=size(circle_2,1)  
    if size(circle_2{i,1},1)<10
        circle_2(i,:)=[];
        i=i-1;
    end
    i=i+1;
end
%% 将circle_1和circle2中拟合参数存在la中
la=[];l=1;
%存入断裂点拟合圆参数
for i=1:1:size(circle_1,1)
    [p,~]=Circle_Fitting(circle_1{i,1} (:,2:3));
    la(l,1)=p(1); %横坐标
    la(l,2)=p(2); %纵坐标
    l=l+1;
end

%存入恢复点拟合圆参数
for i=1:1:size(circle_2,1)
    [p,~]=Circle_Fitting(circle_2{i,1} (:,5:6));
    la(l,1)=p(1); %横坐标
    la(l,2)=p(2); %纵坐标
    l=l+1;
end
%% 参数合并
location=[]; %circle_1存储合并后的数据
i=1;k=1;num=0;
while i<=size(la,1)
    location(i,:)=la(i,:);
    j=i+1;
    while j<=size(la,1)
        center1=la(i,1:2);
        center2=la(j,1:2);
        if pdist([center1;center2])<10 
            location(i,1)=(center1(1)+center2(1))/2;
            location(i,2)=(center1(2)+center2(2))/2;
            la(j,:)=[];
            j=j-1;
            num=1;
        end
        j=j+1;
    end
    i=i+1;
end
if isempty(location) %如果为空
    T=4;
    location_detect=[];LocationError_avg=[];point1=[];point2=[];
    return 
end
%% 计算误差
%真实干扰源位置排序
location_real=sortrows(location_real);
%真实干扰源的KM坐标
Interf_KM=[];
Interf_KM(:,1)=(location_real(:,1)-lon1).*(111*cos(location_real(:,2).*(pi/180)));
Interf_KM(:,2)=(location_real(:,2)-lat1).*111;
Interf_KM=sortrows(Interf_KM); %真实干扰源KM位置
location=sortrows(location); %预测干扰源位置排序
[Interf_KM,location,LocationError_avg] = Error(Interf_KM,location);

%% 干扰源位置存储
location2_detect=location(:,2)/111+lat1;
location1_detect=location(:,1)./(111*cos(location2_detect*pi/180))+lon1;
location_detect=[location1_detect,location2_detect,location(:,3)];

location_real = [location_real,Interf_KM(:,3)]; %干扰源真实位置
%% 函数
function [InterferencePostion,InterferencePostion_detected,LocationError_avg] = Error(InterferencePostion,InterferencePostion_detected)
InterferenceNum = size(InterferencePostion,1);%真实干扰源数量
InterferenceNum_detected = size(InterferencePostion_detected,1);%干扰源检测数量
LocationErrorr=zeros(1,1);
for mm= 1: InterferenceNum_detected
    for n = 1: InterferenceNum
        LocationErrorr(mm,n)=pdist([InterferencePostion_detected(mm,1),InterferencePostion_detected(mm,2);InterferencePostion(n,1),InterferencePostion(n,2)]);
    end
end
if InterferenceNum_detected>=InterferenceNum
   [LocationError_min,I] = min(LocationErrorr,[],1);   %计算每个干扰源定位结果到干扰源位置的最小距离
   InterferencePostion_detected(:,3)=-1; %虚检
   InterferencePostion_detected(I,3)=1; %无虚检
   InterferencePostion(:,3)=LocationError_min'; %定位误差
end
if InterferenceNum_detected<InterferenceNum
   [LocationError_min,I] = min(LocationErrorr,[],2);   %计算每个干扰源定位结果到干扰源位置的最小距离
   InterferencePostion_detected(:,3)=1; %无虚检
   InterferencePostion(:,3)=-1; %漏检
   InterferencePostion(I,3)=LocationError_min';  %无漏检定位误差
end
  LocationError_avg = mean(LocationError_min);   %计算所有干扰源定位结果的平均误差
end

function Z=Gauss2D(X,Y,sigma,a,b)
%XY是网格坐标
%sigma是高斯分布的宽度
%ab是中心点坐标
    Z=0.5*exp(-((X-a).^2+(Y-b).^2)./sigma.^2);
end

function B=Heatsum(A1,A2)
%两个点之间叠加
B=1-(1-A1).*(1-A2);
end
%% 计算误差参数m的均值
function t1=aver_m(circle)
 t1=0;
 for ss=1:size(circle,1)
     t1=t1+circle{ss,3};
 end
 t1=t1/size(circle,1);
end
%% 断裂点集拟合圆，数据存储在circle中
function [circle]=cirlce(Data,s)
circle={};
numm=0;
%每个点集拟合圆
circle={}; %存储拟合圆数据点集
numm=0;
if s==1
    DAta=Data(:,2:3);
    DAta_1=Data(:,2:3);
end
if s==2
    DAta=Data(:,5:6);
    DAta_1=Data(:,5:6);
end

    %第一段拟合圆
    [pp,mm]=Circle_Fitting(DAta_1);
    while mm>0.02 %0.02
        DAta_1=DAta_1(1:end-1,:);
        [pp,mm]=Circle_Fitting(DAta_1);
    end
    len1=size(DAta_1,1);
    numm=numm+1;
    circle{numm,1}=Data(1:len1,1:7);
    circle{numm,2}=pp;
    circle{numm,3}=mm;
    %第二段及以后的拟合圆
    while (len1+2)<size(DAta,1)
        numm=numm+1;
        data_2=DAta(len1+2:end,:);
        [pp,mm]=Circle_Fitting(data_2);
        while mm>0.02 %0.02
            data_2=data_2(1:end-1,:);
            [pp,mm]=Circle_Fitting(data_2);
        end
        len2=size(data_2,1);
        circle{numm,1}=Data(len1+2:len1+len2+1,1:7); %后面的拟合圆
        circle{numm,2}=pp;
        circle{numm,3}=mm;
        len1=len1+len2+1; 
    end
end
%% 三点定圆
function [circleCenter,radius] = ThreePointCircle(x,y,z)
%circleCenter：圆心
%radius半径
%x,y,z为三个点
            A=[x(1)-y(1),x(2)-y(2);z(1)-y(1),z(2)-y(2)];
            B=[x(1)^2-y(1)^2+x(2)^2-y(2)^2;z(1)^2-y(1)^2+z(2)^2-y(2)^2];
            
            ab=A\B;
            aa=ab(1)/2;
            bb=ab(2)/2;
            circleCenter = [aa,bb];

            c2 = (x(1)-aa)^2+(x(2)-bb)^2;
            radius = sqrt(c2);
end
%% 最小二乘拟合圆的函数
function [p,m] = Circle_Fitting( XZ )
N = size(XZ,1);
x = XZ(:,1);
z = XZ(:,2);
sum_X_Raw = 0;
sum_Z_Raw = 0;
sum_XSquare_Raw = 0;
sum_ZSquare_Raw = 0;
sum_XCube_Raw = 0;
sum_ZCube_Raw = 0;
sum_XZZ_Raw = 0;
sum_XZ_Raw = 0;
sum_XXZ_Raw = 0;

for ii=1:N
    sum_X_Raw = sum_X_Raw+x(ii);
    sum_Z_Raw = sum_Z_Raw+z(ii);
    sum_XSquare_Raw = sum_XSquare_Raw+x(ii)*x(ii);
    sum_ZSquare_Raw = sum_ZSquare_Raw+z(ii)*z(ii);
    sum_XCube_Raw = sum_XCube_Raw+x(ii)*x(ii)*x(ii);
    sum_ZCube_Raw = sum_ZCube_Raw+z(ii)*z(ii)*z(ii);
    sum_XZ_Raw = sum_XZ_Raw+x(ii)*z(ii);
    sum_XZZ_Raw = sum_XZZ_Raw+x(ii)*z(ii)*z(ii);
    sum_XXZ_Raw = sum_XXZ_Raw+x(ii)*x(ii)*z(ii);
end
D = N*sum_XZ_Raw-sum_X_Raw*sum_Z_Raw;
C = N*sum_XSquare_Raw-sum_X_Raw*sum_X_Raw;
E = N*sum_XCube_Raw+N*sum_XZZ_Raw-(sum_XSquare_Raw+sum_ZSquare_Raw)*sum_X_Raw;
G = N*sum_ZSquare_Raw-sum_Z_Raw*sum_Z_Raw;
HH = N*sum_ZCube_Raw+N*sum_XXZ_Raw-(sum_XSquare_Raw+sum_ZSquare_Raw)*sum_Z_Raw;

aa = (HH*D-E*G)/(C*G-D*D);
bb = (HH*C-E*D)/(D*D-G*C);
c = -((sum_XSquare_Raw+sum_ZSquare_Raw)+aa*sum_X_Raw+bb*sum_Z_Raw)/N;

p(1) = -0.5*aa;
p(2) = -0.5*bb;
p(3) = 0.5*sqrt(aa*aa+bb*bb-4*c);
% m=max(abs(x.^2+z.^2+x.*a+z.*b+c));
mm=0;
for ii=1:1:size(XZ,1)
    mm=mm+abs(pdist([p(1:2);XZ(ii,1:2)])-p(3));
end
m=mm/size(XZ,1)/p(3);
end
%% 排序,对每一个点集排序
function clusterData1=sequence(clusterData1,s)
%计算中心点
x=[];y=[];
for ii=1:size(clusterData1,1)
    dataa=clusterData1{ii,1};
    if s==1
        x=[x;dataa(:,2)];y=[y;dataa(:,3)];
    end
    if s==2
        x=[x;dataa(:,5)];y=[y;dataa(:,6)];
    end  
end
cx = mean(x);
cy = mean(y);
%排序,对每一个点集排序
for ii=1:size(clusterData1,1)
    dataa=clusterData1{ii,1};
    if s==1
        x=dataa(:,2);y=dataa(:,3);
    end
    if s==2
        x=dataa(:,5);y=dataa(:,6);
    end
    aa = atan2(y - cy, x - cx); %求角度
    [~, order] = sort(aa);%找到正确排序
    for tt=1:size(dataa,2)
        dataa(:,tt)=dataa(order,tt);
    end
    clusterData1{ii,1}=dataa;
end
end


end